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полной матрицы 


В работе предлагается подход, который базируется на неявных многостадийных методах, модифицированных 
таким образом, что неявные стадии становятся параллельными. Сокращение числа обменов достигается за 
счет преобразования исходной матрицы, приводящего функциональный определитель к диагональному 
виду. Получаемые на основе такого подхода расчетные схемы обладают меньшей вычислительной сложностью 
и являются весьма эффективными при решении жестких уравнений. 


Введение 


Данная работа является продолжением исследований, представленных в [1-6], 
и посвящена параллельной численной реализации решения задачи Коши с помощью 
неявных методов, ориентированных на решение жестких уравнений и их систем. 
Интегрирование таких уравнений основано на формировании и решении на каждом 
шаге нелинейной системы алгебраических уравнений размерностью тх 5, где т — 
размерность системы или наивысший порядок уравнения, 5 — число стадий метода. 
Эффективное решение такой системы является главной проблемой при реализации 
неявного стадийного метода [7-10]. Классический подход заключается в использовании 
итерационного метода Ньютона с полным якобианом [7]. При решении жестких 
систем этот подход является непривлекательным из-за высокой трудоемкости реали- 
зации, связанной, прежде всего, с необходимостью многократного переопределения 
величины шага интегрирования на участках быстрого изменения производной [1], [3], [4]. 

Для решения жестких систем обыкновенных дифференциальных уравнений на 
сегодняшний день предлагаются различные типы параллельных методов [8-10]. Большая 
часть их базируется на неявных методах, обычно классических многостадийных методах 
Рунге-Кутты, в которых неявные отношения разрешаются с помощью итерационного 
процесса. В пределах каждой итерации стадийный метод обладает хорошим парал- 
лелизмом, так как вычисление компонентов вектора на итерации распределяется по 5 
процессорам. Однако после каждой итерации процессоры должны обмениваться 
полученными результатами, а это подразумевает частую связь между процессорами. 
Такой мелкозернистый параллелизм особенно не привлекателен при использовании 
в компьютерах с распределенной памятью. 

В работе предлагается альтернативный подход, который также базируется на 
неявных многостадийных методах, модифицированных таким образом, что неявные 
стадии являются уже параллельными, так, что значения в стадийных точках могут быть 
получены независимо друг от друга. То есть обмен значениями процессоры осу- 
ществляют не после каждой итерации, а после получения значения для очередной 
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расчетной точки. Такое радикальное сокращение числа обменов достигается за счет 
использования диагонального приближения исходной матрицы, приводящее исходный 
функциональный определитель к диагональному виду. Получаемые на основе такого 
подхода расчетные схемы обладают меньшей вычислительной сложностью и являются 
весьма эффективными при решении жестких уравнений. 

Цель данной работы состоит в создании расчетных схем для параллельного 
решения жестких уравнений и их систем, приводящих исходный функциональный 
определитель к диагональному виду, что обеспечивает сокращение числа обменов на 
итерациях по стадийным точкам. 


Диагонализация исходной матрицы неявного метода 


При реализации численного решения задачи Коши 


фх 
И» хп) = хо (1) 
с помощью неявных методов 
с| А (2) 
ы 


полная матрица А предварительно приводится к диагональному виду. При этом 
среди множества возможных вариантов приближения Р выбираются такие, что 
спектральный радиус 


вИ-Б”А)=0 (3) 
ИЛИ 
РА-Б'А) > тт. (4) 


Решение задачи, связанной с поиском корней нелинейного уравнения (3), 
эквивалентно задаче 


ае(р"А- АГ) =0, (5) 
количество корней которой определяется стадийностью метода и равно 5: 
аеКрА- М) =(1-А)”. 


Это уравнение преобразовывается к следующему виду 
(=) + У и(-д м = + У СК-А У, (6) 
1=1 1=1 


где С, — комбинаторное соединение числа сочетаний, 


и; - коэффициенты характеристического многочлена РА. 


и; оцениваются как сумма миноров 1-Й степени, симметричных относительно 


главной диагонали РА. В частности, предельные случаи 
И! = пасе( 7! А), и; = ае(Р-'А). 


В (6) левая и правая части эквивалентны, если соответствующие коэффициенты 

и; и С’ совпадают. На этом основании формируется нелинейная система из 5 
неизвестных, такая что 

и: =С', 1=1,2,..5. (7) 
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Из множества возможных решений этой системы, максимальное число которых 


определяется как 2” — 1, выбирается такое 
ам азиая 


на котором выполняются соотношения (3) или (4). Этот вариант Р-! обеспечит 
лучшую сходимость численной реализации. Для такого выбора на каждом полученном 
решении строится матрица вида 


1 
1(2)=2(1-20) (А-Б) (8) 
и исследуется спектр ее собственных значений. По каждому полученному решению 
оценивается максимальное собственное значение тах(р(7(2)) спектра матрицы 
7(2). Среди множества решений для реализации выбирается вариант с минимальным 


значением тах(р(7(2)). Таким образом, подход для выбора матрицы О, рассматри- 
ваемый в этой работе, базируется на минимизации спектрального радиуса матрицы. 


Модификация неявных методов для жестких задач 
на основе диагональной матрицы 


При построении адекватных математических моделей учитывается большое число 
факторов, что неизбежно приводит к явлению жесткости и описывающим его жестким 
системам. Под жесткими уравнениями будем подразумевать такие, для которых опреде- 
ленные неявные методы дают лучший результат, обычно несопоставимый с явными 
методами [7]. Формализация такого прагматического определения для разных типов 
уравнений и систем может носить различный характер [8]. 


[и-1 -1+САЙн [-1+С2Й» нс =Ь 


и 
Й» 


Рисунок 1 — Шаблон 5 — стадийного неявного метода 


В качестве исходного метода решения жесткой системы выбирается 5 — стадий- 
ный неявный метод 


СПЕ Х +В (АХЕ(Х()), 1=12,...,4, 
х: Ех +И,(Ь х1.)Е( СЯ? ) (9) 
с шаблоном, приведенным на рис. 1. Параллельный вариант реализации такого метода 
ХО (ах) + (ржа Ех), }=12,..4, 
хх, +, (ВТ хГ,)Е(Х ) (10) 


с шаблоном, приведенным на рис. 2, модифицируется. 
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[и-1 [..1-С1Йн [,-1-С2Йя 1-Й = 


у 
Й» 


Рисунок 2 — Шаблон 5 — стадийного параллельного метода 


Для этого вводятся дополнительные стадийные точки, которые были получены 
на последней итерации при расчете значения для точки {„_, следующим образом 


х(6) = Р(Х(@) м й = /[,2,..., М, 


п-1, 


Шаблон расчетной схемы приведен на рис. 3. Тогда модифицированный параллель- 
ный полностью неявный метод будет иметь вид 


ХИ -в(РхГ,)Е(Х(” )=ехх_ +1 ((А-Р)х Г, )Е(Х”), }=12,...,4, 
же +й, (6х1, )Е(Х (9), (11) 


где 5 обозначает количество стадий неявного стадийного метода, описанного с 
помощью с, 4, ир (2). 


[„2 24-Й | 1-Й [-1+С2Йн [1+5 Ин-Е 


у 
1 Й, 


Рисунок 3 — Шаблон модифицированного 5 — стадийного параллельного метода 


Вектор Х„ состоит из 5 вложенных стадийных векторов х„/,х„2,..,Х„.‹, Каждый 
являющийся приближением к решению в промежуточных стадиях по времени 
1+сй,, и Е(Х„) состоит из вложенных векторов Л(х„1), /(Хи2),..., Л(Хиз Л. 


[; обозначает единичную матрицу размерности 5 ие = (1..., 1) - единичный вектор 


размерности 5. Р — вектор стадийных значений последней итерации при расчете 
значения для {„_;, М - количество шагов по времени, и 4 — количество итераций, 


обеспечивающее вычисление корректора с заданной точностью. Исходя из того, что 
р является диагональной матрицей, 5 стадийных векторов в Х 0 й могут быть 


вычислены параллельно. 
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Исследование устойчивости неявных методов 
на основе диагональной матрицы 


При исследовании устойчивости блочных разностных методов для жестких 
систем уравнений, так же, как и для классических методов, обычно рассматривают 
модельное уравнение 

хе = Ях, (12) 
МА 

где — произвольное комплексное число. Свойства различных методов анализируют 
на примере модельного уравнения (12). Для того чтобы уравнение (11) действительно 
моделировало исходную систему (1), необходимо рассматривать его при всех таких 
Я, которые являются собственными числами матрицы (8). Кроме того, все корни 
характеристического уравнения (8) не должны превосходить по модулю единицу. 
При |2] — 0 накапливаются нежесткие компоненты ошибки 


7(2)>2(А-О), 
при |2] —> © — жесткие компоненты ошибки 
7(2)>1-Р!А. 


Собственные числа матрицы /-— Р-'А при этом должны соответствовать усло- 
виям (3-4). Плохая обусловленность собственных чисел — одна из причин, усложняющая 
минимизацию спектрального радиуса матрицы и увеличивающая время получения 


решения. В работе предлагается не только минимизация (1 -Р-'А) для неявных ме- 


у ЕТ 
тодов с числом стадий 5 > 2, но и поиск точного решения р(1-Р А) =0. Использо- 


вание Майетайса для символьных вычислений значительно упрощает эту процедуру. 
При этом можно получать как точные решения нелинейных уравнений, что возможно при 
5 < 4 [8], так и численные решения с заданной степенью точности. 


Численная реализация построения диагональной матрицы 


Рассмотрим несколько вариантов построения диагональных матриц Ш) для 
многостадийных методов. В качестве исходных для диагонализации выбираются 
известные полностью неявные стадийные методы. Для каждого из выбранных 
исходных методов определяется лучший вариант диагонализации. 

Выберем в качестве варианта диагонализации неявный трехстадийный метод 
Лобатто ШС порядка 4 [7]. Из исходной матрицы А 


7015 0 
176 Э/1а. =41/12 (13) 
1/6 2/3 1/6 


формируется матрица Р-”А с искомыми элементами О”! = фае(4 »„@2, 43) 


а а 4, 
63 б 
а, 34а, а), 
6 12 12 
о 
63 бб 
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Приравнивание коэффициентов полинома в соответствии с (6) — (7) приведет к 
следующей нелинейной системе 


4, 54, 4 
= И 
о ша 
4142 ‚ 4243 
+ =З, 
И? 8 8 
муар“) “2-1. 


Решением такой системы будут являться 4 действительных корня. Лучшему ва- 
рианту диагонализации соответствует минимум из всех тах(р(2(2))) = 0.5096384 . 


ши т 
0.3 + 


0.2 . 
0.1: 


0.2 . 


0.1: 


0.1 0.2 0.3 0.4 0.5 0.6 0 02 03 04 ор 


0.11 
0.2 + 
0.3 + 


0.1: 


0.2 . 


а) 6) 
Рисунок 4 — Годографы функции 2(2) для метода Лобатто ШС порядка 4 


Оценки параллелизма разработанных методов 


Реализация предлагаемых алгоритмов, основанных на диагональных преобра- 
зованиях, ориентирована на использование многопроцессорных вычислительных систем 
с линейкой процессорных элементов. Набор процессоров известен до начала вычислений 
и не меняется в процессе счета, при этом каждый процессорный элемент может 
выполнить любую арифметическую операцию за один такт, временные затраты, связан- 
ные с обращением к запоминающему устройству, отсутствуют. В качестве примера рас- 
смотрим неявный трехстадийный метод Лобатто ШС порядка 4 с исходной матрицей (13), 
для которой был получен лучший вариант диагонализации 


0.6537 0 0 
0 0.184059 0 
0 0 0.3462998 


Характеристики параллелизма, ускорение и эффективность, исследовались для 
систем с изменяющимися трудоемкостями вычисления правых частей №, прини- 
мающими значения #={ 10, 50, 100, 500, 5000}. Реализация двухстадийного неявного 
метода с нижней треугольной матрицей на 5ПМО структуре с числом процессорных 
элементов, совпадающих с размерностью системы т, дает следующие показатели 
ускорения и эффективности. 

Сравнение полученных показателей по группам позволяет утверждать, что по- 
казатели неявных методов с диагональными матрицами значительно превосходят 
соответствующие характеристики для методов полными матрицами, при этом по- 
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казатели тем лучше, чем выше трудоемкости вычисления правых частей, что 
объясняется сокращением отношения общего времени счета к времени реализации 
обменов между процессорными элементами. 


20 40 60 80 10 


Рисунок 5 — Характеристики параллелизма 3-стадийного метода Лобатто ШС 
с диагональной матрицей 


Выводы 


Работа направлена на сокращение числа обменов при параллельной численной 
реализации решения задачи Коши с помощью неявных методов, ориентированных на 
решение жестких уравнений и их систем. Предлагаемый подход базируется на 
модификации неявных многостадийных методов, обеспечивающей параллельное 
получение значений в стадийных точках, при этом обмен значениями процессоры 
осуществляют не после каждой итерации, а после получения значения для очередной 
расчетной точки. Такое радикальное сокращение числа обменов достигается за счет 
использования диагонального приближения исходной матрицы. Эффективность решения 
такой системы обеспечивается диагонализацией исходной матрицы неявного метода с 
обеспечением максимальной скорости сходимости при параллельной реализации. 

Для классических неявных методов, которые используются для решения жестких 
уравнений и их систем, диагональные вхождения 2) сложны и требуют дальнейших 
модификаций, вовлекающих в расчеты сложную арифметику. Но окончательно постро- 
енные таким образом итерационные методы устойчивы, характеризуются высокой 
скоростью сходимости, обладают естественным параллелизмом, что обеспечивает их 
эффективную реализацию в параллельных вычислительных системах. 

Рассмотрены варианты диагонализации известных стадийных методов. Приведена 
методика выбора лучшего варианта диагонального приближения, основывающаяся на 
минимизации спектрального радиуса матрицы. 
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оенчм 


О.А. Дмитруева 

Паралельне моделювання жорстких систем на основ! дагоналзацй повно! матрищ 

У робот! пропонуеться шджд, що базуеться на неявних багатостаминих методах, модификованих таким 
чином, що неявн! стадй стають паралельними. Скорочення числа обман в досягаеться за рахунок перетворення 
видно! матриц що приводить функцюнальний визначник до дагонального виду. Одержуван! на основ! 
такого шдходу розрахунков! схеми мають меншу обчислювальну складысть 1 е досить ефективними при 
розв’язанн! жорстких равнянь. 


О.А. Ртитчеуа 

РагаП@а МодеНпо оЁ ЗИЁ Зузетб оп Фе Ва$5 оГа П1агопа|7айоп оГа Еи| Маблх 

ТЬЫ5$ \огК 51296545 Фе арргоасВ, \сЬ 1$ Базе4 оп е парИсй ши@яер тефо4$, по 1Йе4 11 засВ а \ау 
Фаё парПси 5асез Бесоште рагаПе|. Кедасйоп оЁ Фе питбег оЁ ехсБапое$ 1$ геасВе4 дие ю тапзРогтаноп 
оЁ пе шша| тафлх, [еа4те а Епсйопа| де4егилтап: 0 а Ф1асопа| Кта. Оез1юп зсВетез саште оп фе Ба$1$ 
ОЕ зисВ арргоасЬ Бауе зта|ег сотрийпе сотр]ехиИу ап4 аге гафег еЙесйуе ш зо!уше 5ИЁР едиайопз. 


Статья поступила в редакцию 23.08.2011. 
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